/*
 * pegasos.cpp
 *
 *  Created on: Nov 7, 2013
 *      Author: jieshen
 */

#include "cblas.h"
#include <cstdio>
#include <cstdlib>
#include <utility>
#include <vector>
#include <iostream>

using namespace std;

/*
 * w: weight
 * dim: dimension
 * lambda and K, refer to the pegasos paper
 * F: samples: pair<feature, label>
 */

void update_model_with_pegasos(float* w, const int dim,
                               const float lambda /*=0.001*/,
                               const size_t K /*=1*/,
                               const vector<pair<vector<float>*, int> >& F)
{
  int iter = 0;

  const size_t sz = sizeof(float) * dim;

  vector<size_t> idx(F.size(), 0);
  for (size_t i = 0; i < idx.size(); ++i)
    idx[i] = i;

  const float w_norm2 = cblas_snrm2(dim, w, 1);
  const float fact = std::min(1.0, 1.0 / (std::sqrt(lambda) * w_norm2));
  cblas_sscal(dim, fact, w, 1);

  const int max_iter = 100 / lambda;

  float* sum_hard_samples = (float*) malloc(sz);
  memset(sum_hard_samples, 0, sz);

  float* w_semi = (float*) malloc(sz);
  memset(w_semi, 0, sz);

  cerr << "\tPegasos update iteration: ";
  fflush(stderr);

  clock_t begin = clock();

  while (iter++ < max_iter)
  {
    if (iter % (max_iter / 100) == 0)
    {
      cerr << ".";
      fflush(stderr);
    }

    const float eta = 1.0 / (lambda * (iter + 1));

    std::random_shuffle(idx.begin(), idx.end());

    memset(sum_hard_samples, 0, sz);

    for (size_t i = 0; i < K; ++i)
    {
      const pair<vector<float>*, int>& cur_sample = F[idx[i]];
      if ((size_t) dim != cur_sample.first->size())
      {
        cerr << dim << " " << cur_sample.first->size() << endl;
        exit(0);
      }
      const float pred = cblas_sdot(dim, w, 1, cur_sample.first->data(), 1);
      const float hinge_loss = std::max(1.0 - cur_sample.second * pred, 0.0);

      if (hinge_loss > 0)
        cblas_saxpy(dim, cur_sample.second, cur_sample.first->data(), 1,
                    sum_hard_samples, 1);
    }

    memset(w_semi, 0, sz);
    cblas_saxpy(dim, eta / K, sum_hard_samples, 1, w_semi, 1);
    cblas_saxpy(dim, 1.0 - eta * lambda, w, 1, w_semi, 1);

    const float w_semi_norm2 = cblas_snrm2(dim, w_semi, 1);
    const float factor = std::min(1.0,
                                  1.0 / (std::sqrt(lambda) * w_semi_norm2));

    memset(w, 0, sz);
    cblas_saxpy(dim, factor, w_semi, 1, w, 1);

    //cerr << "iter=" << iter << "done" << endl;
  }

  free(sum_hard_samples);
  free(w_semi);

  clock_t end = clock();

  cerr << (end - begin) / CLOCKS_PER_SEC << " seconds" << endl;
}

